{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e174828a",
   "metadata": {},
   "source": [
    "## Changing Temperature and Pressure\n",
    "\n",
    "Sometimes you want the temperature and/or pressure of a simulation to vary with time.  To see how to do this, let's begin by loading a PDB file and creating a Simulation for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3592a17b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from openmm.app import *\n",
    "from openmm import *\n",
    "from openmm.unit import *\n",
    "\n",
    "pdb = PDBFile('villin.pdb')\n",
    "forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')\n",
    "system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,\n",
    "                                 nonbondedCutoff=1*nanometer, constraints=HBonds)\n",
    "integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)\n",
    "simulation = Simulation(pdb.topology, system, integrator)\n",
    "simulation.context.setPositions(pdb.positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4819f91e",
   "metadata": {},
   "source": [
    "When we created the LangevinMiddleIntegrator, we specified the temperature as `300*kelvin`.  We are free to change that at any time.  Let's run a short simulation, steadily decreasing the temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fbe9edaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(10):\n",
    "    integrator.setTemperature((300-30*i)*kelvin)\n",
    "    integrator.step(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c836796",
   "metadata": {},
   "source": [
    "Things are a little more complicated when we have a barostat.  The Monte Carlo acceptance criterion for the barostat depends on temperature, so we need to tell both the integrator and the barostat about the new temperature.\n",
    "\n",
    "Let's add a barostat to the System.  Note the call to `reinitialize()`.  Without that, the existing Simulation will not see any changes we make to the System."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0dd62ed9",
   "metadata": {},
   "outputs": [],
   "source": [
    "system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))\n",
    "simulation.context.reinitialize(preserveState=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3a370dc",
   "metadata": {},
   "source": [
    "Because the barostat is part of the System, it stores all its state information in the Context.  We can change that information with `setParameter()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bab6dd29",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(10):\n",
    "    temperature = (300-30*i)*kelvin\n",
    "    simulation.context.setParameter(MonteCarloBarostat.Temperature(), temperature)\n",
    "    integrator.setTemperature(temperature)\n",
    "    integrator.step(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f3d41ca",
   "metadata": {},
   "source": [
    "Setting the pressure works exactly the same way, just with a different parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "978be843",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.context.setParameter(MonteCarloBarostat.Pressure(), 0.5*bar)"
   ]
  }
 ],
 "metadata": {
  "tags": ["barostat", "thermostat", "application layer"],
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}